Export fasta

  cultures_Pelago_18S <-  cultures_Pelago %>%
    filter(!is.na(accession_18S)) %>% 
    mutate(sequence_name = str_c(accession_18S, str_c("RCC", rcc_id), strain_name,species, sep="|" ))
  
  seq_out <- DNAStringSet(cultures_Pelago_18S$sequence_18S)  # Store the sequence in a  DNAString
  names(seq_out) <- cultures_Pelago_18S$sequence_name
  
  writeXStringSet(seq_out, "C:/Users/vaulot/Desktop/RCC_Pelago_18S.fasta", compress=FALSE) 

Update and complement data

• Add combination of samples and species that have 0 reads
• Add class specific information

update_asv_set <- function(one_asv_set) {
  
# one_asv_set <- asv_set[["V4"]]
  
one_asv_set$samples <- one_asv_set$samples %>% 
    filter(!is.na(file_code) ) %>% 
    filter(file_code !="") %>% 
    filter(!is.na(DNA_RNA))  # Some Tara files do not correspond to any samples...
  

# Find samples which have no pelagophytes  
samples_no_pelago <- one_asv_set$samples %>% 
    filter(!(file_code %in% one_asv_set$df$file_code))  

# Extract just file_code, species and n_reads
samples_species_counts <- one_asv_set$df %>% 
  select(file_code, species, n_reads)
taxo <- one_asv_set$df %>% 
  select(species, order, family, genus, species) %>% 
  distinct()


# Create a full table with all combinations of samples and species

samples_species <- bind_rows(one_asv_set$df, samples_no_pelago) %>% 
  expand(file_code, species) %>% 
  filter(!is.na(species)) %>% 
  left_join(one_asv_set$samples) %>% 
  left_join(taxo) %>% 
  left_join(samples_species_counts) %>% 
  mutate(n_reads = replace_na(n_reads, 0))

one_asv_set$df <- samples_species

  
# Normalize by total number of Pelagophyceae reads
  asv_class <- one_asv_set$df %>% 
    group_by(file_code) %>% 
    summarise(reads_class = sum(n_reads)) 
  
  one_asv_set$df <- left_join(one_asv_set$df, asv_class) %>% 
    mutate (pct_class = n_reads/reads_class*100,
            pct_total = n_reads/reads_total*100,
            pct_pelago_total = reads_class/reads_total*100)
  
  return(one_asv_set)
  }

Statistics

summarise_asv_set <- function(one_asv_set, directory) {
  
# one_asv_set <- filter(asv_set[["V4"]]$df, dataset_id != 36)
# directory = "../metabarcodes/V4/"
  
  summary_species <- phyloseq_long_treemap(one_asv_set$df, genus, species, "Species" )

  
  summary_pct <- one_asv_set$df %>% 
    group_by(species) %>% 
    summarise(pct_class = mean(pct_class, na.rm = TRUE),
              pct_total = mean(pct_total, na.rm = TRUE))
  
  summary_species <- left_join(summary_species$df, summary_pct)
  
  openxlsx::write.xlsx(summary_species, str_c(directory, "summary_species.xlsx"))
  

  
  summary_substrate <- one_asv_set$samples %>% 
    count(substrate)
  
  summary_depth <- one_asv_set$samples %>% 
    count(depth_level) 
  
  summary_fraction <- one_asv_set$samples %>% 
    mutate(fraction_label = str_c(format(fraction_min, nsmall=2, scientific = FALSE), 
                                  format(fraction_max, nsmall=2, scientific = FALSE), sep="_") ) %>% 
    count(fraction_label) 
  

  return(list(summary_species, summary_substrate, summary_depth, summary_fraction))
   
    
  }

Maps

map_asv_set <- function(one_asv_set, directory) {
  
# one_asv_set <- filter(asv_set[["V4"]]$df, dataset_id != 36)
# directory = "../metabarcodes/V4/"
  
  world <- map_data("world")

  # g <- ggplot() + 
  #   geom_polygon(data = world, aes(x=long, y = lat, group = group), fill="grey") + 
  #   coord_fixed(1.3) +
  #   geom_point(data=one_asv_set, aes(x=longitude, y=latitude, size=pct_class), fill="blue", shape=21) + 
  #   # scale_size_area(trans="log10") +
  #   ggtitle("Percent of Pelagophyceae") +
  #   facet_wrap(vars(species))
  # print(g)

  df <- one_asv_set  %>% 
    select(longitude, latitude, pct_pelago_total) %>% 
    rename(z=pct_pelago_total) %>% 
    distinct()
  gg <- map_distribution(df, 
                         map_title = "Pelagophyceae - % of Total reads",
                         z_limits = c(0, 10),
                         z_breaks = 1*c(1, 2, 4, 6, 8, 10),
                         z_min = 0.01)
  print(gg)
  ggsave(plot= gg , filename=str_c(directory,"map_Pelago_pct_total.pdf"),
     width = 15 , height = 12, scale=1.80, 
     units="cm", useDingbats=FALSE)  
  

  for (one_species in unique(one_asv_set$species)) {
    df <- one_asv_set %>% 
      filter(species == one_species) %>% 
      select(longitude, latitude, pct_class) %>% 
      rename(z=pct_class)
    
    gg <- map_distribution(df, 
                           map_title = str_c(one_species," - % of Pelagophyceae"),
                           z_min = 0.01)
    print(gg)
    ggsave(plot= gg , filename=str_c(directory,"maps_pct_pelago/map_",one_species,"_pct_pelago.pdf"),
       width = 15 , height = 12, scale=1.80, 
       units="cm", useDingbats=FALSE)
  }
  
  # for (one_species in unique(one_asv_set$species)) {
  #   df <- one_asv_set %>% 
  #     filter(species == one_species) %>% 
  #     select(longitude, latitude, pct_total) %>% 
  #     rename(z=pct_total)
  #   gg <- map_distribution(df, map_title = str_c(one_species," - % of Total reads"),
  #                          z_limits = c(0, 1),
  #                          z_breaks = 0.1*c(1, 2, 4, 6, 8, 10))
  #   print(gg)
  #   ggsave(plot= gg , filename=str_c(directory,"maps_pct_total/map_",one_species,"_pct_total.pdf"),
  #      width = 15 , height = 12, scale=1.80, 
  #      units="cm", useDingbats=FALSE)
  # }


  return(TRUE)
   
    
  }

Time series

ts_asv_set <- function(one_asv_set, directory) {
  
  
  # one_asv_set <- filter(asv_set[["V4"]]$df, dataset_id == 36)

  df <- one_asv_set  %>% 
    mutate(replicate = replace_na(replicate, 1)) %>%
    filter(replicate == 1) %>% 
    select(date, species, date, fraction_name, pct_class) %>% 
    mutate(date = lubridate::as_date(date))
  
  
 gg <- ggplot(df, aes(x = date, y =  pct_class, fill = species)) + 
               geom_col(width=30) + 
               theme_bw() + 
               scale_fill_viridis_d() + 
               xlab("") + ylab("Proportion") + 
               facet_grid(rows = vars(fraction_name) ) +
               scale_x_date()
  
  print(gg)
  ggsave(plot= gg , filename=str_c(directory,"Blanes_time_series.pdf"),
     width = 15 , height = 12, scale=1.80,
     units="cm", useDingbats=FALSE)
  

  return(TRUE)
   
    
  }

V4 - OSD 2014-2015, Malaspina, Blanes, Biomarks

  asv_set <- list()
  one_region <- "V4"
  
    asv_set[[one_region]] <- readRDS( str_c("../metabarcodes/", one_region, "/asv_set.rda"))
    
    asv_set[[one_region]] <- update_asv_set(asv_set[[one_region]])
    
    summaries <- summarise_asv_set(asv_set[[one_region]], 
                                   directory = str_c("../metabarcodes/", one_region, "/"))

    purrr::map(summaries,knitr::kable)

[[1]]


genus                      species                         n_reads    pct_class   pct_total
-------------------------  -----------------------------  --------  -----------  ----------
Ankylochrysis              Ankylochrysis_lutea                1124    0.3534281   0.0016651
Ankylochrysis              Ankylochrysis_sp.                   115    0.4789488   0.0002084
Aureococcus                Aureococcus_anophagefferens       32697    7.9762765   0.0310217
Aureoumbra                 Aureoumbra_lagunensis                50    0.1202055   0.0001139
Chrysocystis               Chrysocystis_fragilis              2799    1.4307706   0.0020358
Chrysoreinhardia           Chrysoreinhardia_sp.                  2    0.0000070   0.0000020
Pelagococcus               Pelagococcus_sp.                   6110    3.4071300   0.0075800
Pelagococcus               Pelagococcus_subviridis               1    0.0003316   0.0000011
Pelagomonadaceae_clade_A   Pelagomonadaceae_clade_A_sp.      23710    9.2576480   0.0249426
Pelagomonadaceae_X         Pelagomonadaceae_X_sp.              351    0.1036844   0.0005932
Pelagomonas                Pelagomonas_calceolata           720265   40.5003490   0.6885337
Pelagophyceae_Svalbard     Pelagophyceae_Svalbard_sp.          498    0.2980440   0.0011230
Pelagophyceae_XXX          Pelagophyceae_XXX_sp.             50968   11.2727535   0.0481332
Pulvinaria                 Pulvinaria_sp.                        4    0.0001133   0.0000009
Sarcinochrysidaceae_X      Sarcinochrysidaceae_X_sp.           800    1.5698927   0.0015097
Sarcinochrysis             Sarcinochrysis_sp.                12502    7.6793874   0.0171452

[[2]]


substrate       n
----------  -----
sediment       24
water        1000

[[3]]


depth_level       n
-------------  ----
anoxic            7
bathypelagic     63
bottom           24
DCM              45
euphotic         37
mesopelagic      60
surface         788

[[4]]


fraction_label      n
---------------  ----
NA_     NA         24
0.00_   0.20        9
0.20_   3.00      120
0.22_     NA      293
0.22_   3.00      303
0.80_   3.00       38
3.00_  20.00      205
20.00_2000.00      32

# Maps  of surface samples  
    map_asv_set(filter(asv_set[[one_region]]$df,
                       dataset_id != 36,
                       depth_level=="surface"),
                directory = str_c("../metabarcodes/", one_region, "/"))

[1] TRUE

# Blanes time series
    ts_asv_set(filter(asv_set[[one_region]]$df,
                       dataset_id == 36),
                directory = str_c("../metabarcodes/", one_region, "/"))

[1] TRUE

V9 - Tara Oceans

  asv_set <- list()
  one_region <- "V9" 
    asv_set[[one_region]] <- readRDS( str_c("../metabarcodes/", one_region, "/asv_set.rda"))
    
    asv_set[[one_region]] <- update_asv_set(asv_set[[one_region]])
    
    summarise_asv_set(asv_set[[one_region]], 
                directory = str_c("../metabarcodes/", one_region, "/"))

[[1]]
# A tibble: 6 x 5
# Groups:   genus [6]
  genus                   species                    n_reads pct_class pct_total
  <chr>                   <chr>                        <dbl>     <dbl>     <dbl>
1 Aureococcus             Aureococcus_anophageffere~  887561    16.6   0.0800   
2 Pelagococcus            Pelagococcus_sp.             25784     3.67  0.00374  
3 Pelagomonadaceae_clade~ Pelagomonadaceae_clade_A_~   16250     4.24  0.00157  
4 Pelagomonadaceae_X      Pelagomonadaceae_X_sp.       18632     2.77  0.00158  
5 Pelagomonas             Pelagomonas_calceolata     3210461    70.9   0.311    
6 Sarcinochrysidaceae_X   Sarcinochrysidaceae_X_sp.      473     0.222 0.0000379

[[2]]
# A tibble: 1 x 2
  substrate     n
  <chr>     <int>
1 water       884

[[3]]
# A tibble: 3 x 2
  depth_level     n
  <chr>       <int>
1 euphotic      306
2 mesopelagic    70
3 surface       508

[[4]]
# A tibble: 12 x 2
   fraction_label       n
   <chr>            <int>
 1 "  0.22_   3.00"     1
 2 "  0.80_     NA"   112
 3 "  0.80_   3.00"    29
 4 "  0.80_   5.00"   175
 5 "  0.80_  20.00"     8
 6 "  0.80_ 200.00"     2
 7 "  3.00_     NA"    35
 8 "  3.00_  20.00"     1
 9 "  5.00_  20.00"   165
10 " 20.00_ 180.00"   168
11 " 20.00_ 200.00"     8
12 "180.00_2000.00"   180

    map_asv_set(filter(asv_set[[one_region]]$df,
                       dataset_id != 36,
                       depth_level=="surface"),
                directory = str_c("../metabarcodes/", one_region, "/"))

[1] TRUE